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ABSTRACT 


The  behavior  of  Fe-Cr  alloys  under  irradiation  is  in  part  controlled  by  the 
characteristics  of  point  defects  generated  by  high  energy  collision.  Radiation  enhanced 
diffusion  and  radiation  induced  precipitation  are  among  the  mechanisms  that  lead  to 
changes  in  the  microstructure  under  irradiation,  and  are  thus  controlling  effects  such  as 
swelling  and  a’  precipitation.  Point  defects  in  Fe-Cr  alloys  are  diverse  in  nature  due  to 
their  interaction  with  a  variety  of  local  solute  configurations.  Ab  initio  results  indicate 
that  the  magnetic  structure  of  the  alloy  is  critical  in  determining  its  energetics.  The  ability 
to  model  these  properties  with  classic  potentials  is  still  to  be  proven.  In  this  work  a 
detailed  comparison  between  ab  initio  and  classic  values  of  a  variety  of  point  defects 
configurations  is  perfonned,  testing  in  this  way  the  extent  to  which  classic  potentials  can 
be  reliably  used  for  radiation  damage  studies,  and  evaluating  the  dependence  of  point 
defect  formation  energies  on  Cr  concentration. 
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I.  INTRODUCTION 


The  Fe-Cr  alloys  are  considered  for  use  as  structural  materials  for  Generation  IV 
advanced  nuclear  designs.  High  doses  of  neutron  irradiation  cause  vacancies  and 
interstitials  to  form  in  the  steel  matrix.  The  most  accurate  approach  for  evaluating 
formation  energies  of  point  defects  is  numerical  integration  based  on  the  first  principles 
of  quantum  mechanics.  However,  due  to  the  computational  requirements  posed  by  this 
method,  it  is  prohibitive  to  use  it  for  large  scale  simulations.  It  is  therefore  necessary  to 
provide  an  approximation  of  ab  initio  results  using  Molecular  Dynamics  and  Monte  Carlo 
methodologies  coupled  with  an  adequately  descriptive  semi-empirical  many  body 
potential  for  alloys.  The  potential  evaluated  in  this  work  is  based  on  combining  recent  ab 
initio  results  with  thennodynamic  properties  to  describe  formation  energy  as  a  function  of 
local  composition.  This  thesis  seeks  to  answer  the  following  questions: 

1.  Determine  the  extent  to  which  this  potential  adequately  approximates 
results  obtained  by  ab  initio  calculations  in  pure  elements. 

2.  Determine  how  formation  energies  of  point  defects  vary  as  a  function 
of  Cr  concentration;  determine  their  mutual  relationships  and 
stability  in  all  possible  configurations. 

In  order  to  answer  these  questions,  simulations  were  conducted  using  the 
Lawrence  Livermore  National  Laboratory’s  massively  parallel  super-computing 
hardware  and  a  suite  of  publicly  available  and  custom  designed  materials  science 
software  solutions. 

The  calculations  of  this  thesis  show  that  reasonable  predictions  can  be  obtained 
using  the  potential  described  in  this  work.  Formation  energies  of  point  defects  in  pure 
elements  and  as  a  function  of  Cr  concentration  and  heretofore  unexpected  behavior  of 
formation  energies  versus  Cr  concentration  have  been  discovered. 
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II.  BACKGROUND 


A.  HISTORY  OF  NUCLEAR  POWER 

The  use  of  nuclear  power  for  energy  generation  traces  its  origins  to  the  volatile 
times  preceding  World  War  II.  The  theoretical  basis  for  fission  was  laid  by  Otto  Hahn 
and  Fritz  Strassman  in  the  German  periodical  Naturwissenschaften  in  January  1939, 
where  they  reported  that  an  isotope  of  barium  was  produced  by  neutron  bombardment  of 
uranium.  A  knowledge  network  formed  connecting  Otto  Robert  Frisch  and  Lise  Meitner, 
German  scientists  escaping  Hitler’s  rule  to  Denmark,  with  Niels  Bohr  who  was  on  his 
way  to  visit  Albert  Einstein  at  Princeton.  Once  the  information  got  to  Princeton  which 
had  become  a  Mecca  for  theoretical  physics  when  Dr.  Einstein  joined  its  Institute  for 
Advanced  Studies,  it  set  off  a  “chain  reaction”  in  the  scientific  community.  By  the  end  of 
the  year  over  a  hundred  papers  on  the  subject  of  fission  had  been  published  in  Physical 
Review  and  the  concept  of  a  sustainable  chain  reaction  was  grounded  in  theory  [1], 

The  engineering  required  to  sustain  and  control  a  chain  reaction  was,  at  the  time, 
far  from  trivial.  Nonetheless,  the  first  artificial  nuclear  reactor,  Chicago  Pile  1  (CP-1), 
was  built  and  reached  criticality  on  December  2nd,  1942  at  the  underground  racquetball 
courts  at  the  University  of  Chicago  [2]. 


Figure  1.  Drawing  of  Chicago  Pile  1  [3] 


At  the  time  there  was  no  possibility  of  enriching  the  nuclear  fuel,  so  natural 
uranium  had  to  be  used,  which  contains  only  about  0.7  %  of  the  highly  fissionable 
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uranium-235  isotope,  with  the  remainder  consisting  of  the  more  stable  uranium-238.  It 
was  therefore  necessary  to  use  large  quantities  of  highly  purified  graphite  as  neutron 
moderation  material.  The  moderator  was  necessary  to  slow  down  the  high  energy 
neutrons  produced  by  fission  to  an  energy  range  where  they  would  be  much  more  easily 
absorbed  by  U-235  around  .025  eV,  roughly  the  same  temperature  as  that  of  the 
surrounding  material;  hence  the  name  “thermal.”  [4]  The  fission  cross  section  of  U-235 
as  a  function  of  neutron  energy  is  shown  in  Figure  2. 


NEUTRON  CROSS-SECTIONS  FOR  FISSION  OF  URANIUM  AND  PLUTONIUM 


1  bam  =  10-“m2  1  MaV  =  t#  x 

Figure  2.  Neutron  Cross-Sections  for  Fission  of  Uranium  and  Plutonium  [5] 

As  we  can  see  from  Figure  2,  neutrons  produced  by  fission  will  fall  within  the 
energy  range  of  .1  MeV  to  5  MeV,  but  uranium-235  has  the  highest  fission  cross  section 
at  energies  well  below  that  value.  The  “moderator”  is  therefore  used  to  slow  down  the 
fission-produced  neutrons  and  increase  their  probability  of  causing  fissions.  The 
“moderator”  is  usually  a  highly  purified  substance,  generally  deuterium  (in  the  fonn  of 
heavy  water),  beryllium,  or  graphite  for  natural  uranium  reactors.  These  materials  do  not 
absorb  neutrons  easily  (they  have  a  low  absorption  cross  section),  but  their  nuclei  are 
light  enough  that  neutrons  collide  inelastic  ally,  transferring  their  kinetic  energy  to  the 
moderator  and  therefore  slowing  down  in  the  energy  spectrum  towards  the  energies 
which  make  them  more  likely  to  cause  further  fissions  [4].  A  diagram  showing  a  typical 
moderated  chain  reaction  is  shown  in  Figure  3: 
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Once  it  became  possible  to  enrich  uranium,  therefore  increasing  the  percentage  of 
uranium  235  in  the  fuel,  it  became  possible  to  utilize  light  water  for  both  moderation  and 
cooling  of  large  scale  engineering  systems.  This  advance  generally  heralded  the  arrival  of 
the  next  generation  of  nuclear  reactors.  This  generation  of  reactors  was,  for  the  most  part, 
typified  by  pressurized  water  reactors  (PWR)  and  boiling  water  reactors  (BWR)  in  the 
West  and  the  VVER/RBMK  reactors  in  the  former  Soviet  Union  and  satellite  states. 
Canada,  due  to  the  abundance  of  natural  uranium  ore,  opted  for  a  heavy  water  reactor 
design  (CANDU)  which  used  natural  rather  than  enriched  uranium.  All  of  these  systems 
are  considered  to  be  second  generation  nuclear  reactors,  two  common  threads  of  which 
are  the  unique  design  and  features  of  each  reactor  unit  and  their  thermal-neutron  based 
operating  cycle  [6]. 

B.  FUTURE  DIRECTIONS  IN  NUCLEAR  POWER 

The  newly  emerging,  third  generation  of  reactor  designs  are  equipped  with 
advanced  features  such  as  safety  systems  incorporating  passive  energy  dissipation  or 
natural  processes,  simplifying  their  design  and  allowing  them  to  cope  with  malfunctions 
without  the  need  for  operator  action.  Even  though  these  designs  are  still  confined  to  the 
thennal  end  of  the  fission  cross  section,  third  generation  designs  show  some  marked 
improvements  when  compared  to  their  second  generation  cousins.  Some  of  those  are 
listed  below: 

•  a  standardized  design  for  each  type  to  expedite  licensing,  reduce  capital  cost  and 
reduce  construction  time, 

•  a  simpler  and  more  rugged  design,  making  them  easier  to  operate  and  less 
vulnerable  to  operational  upsets, 

•  higher  availability  and  longer  operating  life  -  typically  60  years, 

•  reduced  possibility  of  core  melt  accidents, 

•  minimal  effect  on  the  environment, 

•  higher  bum-up  to  reduce  fuel  use  and  the  amount  of  waste  [7]. 


6 


After  the  Three  Mile  Island  accident,  United  States  imposed  a  self-imposed 
moratorium  on  building  new  nuclear  reactors.  While  third  generation  reactors  have  been 
built  in  Japan,  Taiwan  and  Europe,  this  technology  stagnated  in  the  United  States. 

In  2000,  however,  eleven  world  nations  formed  The  Generation  IV  International 
Forum  (GIF),  a  research  and  development  consortium  tasked  with  leading  the  way 
toward  innovative  nuclear  energy  systems.  Based  on  eight  far-ranging  technology  goals, 
Generation  IV  nuclear  energy  systems  are  aimed  at  achieving  nuclear  energy’s  potential 
worldwide.  The  objective  is  a  new  generation  of  nuclear  energy  systems  that: 

•  advance  nuclear  safety; 

•  address  nuclear  nonproliferation  and  physical  protection  issues; 

•  are  competitively  priced,  and 

•  minimize  waste  and  optimize  natural  resource  utilization. 

Among  the  six  most  promising  designs  for  Generation  IV,  three  are  thermal 
neutron  spectrum  systems  (the  Very-High-Temperature  Reactor  (VHTR)  and 
Supercritical- Water-Cooled  Reactor  (SCWR))  with  coolants  and  temperatures  that  enable 
hydrogen  or  electricity  production  with  high  efficiency,  and  three  are  fast  neutron 
spectrum  systems  (the  Gas-Cooled  (GFR),  the  Fead-Cooled  (FFR),  and  the  Sodium- 
Cooled  (SFR)  fast  reactors)  that  will  enable  more  effective  management  of  actinides 
through  recycling  of  most  components  in  the  discharged  fuel.  [8] 

A  significant  change  in  Generation  IV  relative  to  previous  reactor  generations  is 
the  planned  commercialization  of  fast  spectrum  systems.  A  fast  spectrum  neutron  reactor 
is  a  reactor  in  which  the  chain  reaction  is  sustained  by  fast  neutrons  without  significant 
thennalization.  Such  a  reactor  needs  no  neutron  moderator,  but  must  generally  use  fuel 
that  is  relatively  rich  in  fissile  material  when  compared  to  that  required  for  a  thermal 
reactor.  A  diagram  outlining  the  differences  between  fast  and  thermal  spectrum  processes 
is  shown  in  Figure  4. 
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Figure  4.  Overview  of  Differences  Between  Fast  and  Thermal  Absorption  Processes 


Because  parasitic  absorption  in  the  moderator  can  result  in  a  major  loss  of 
neutrons  in  a  thermal  reactor,  a  fast  reactor  has  an  inherently  superior  neutron  economy; 
that  is,  there  are  excess  neutrons  not  required  to  sustain  the  chain  reaction.  These 
neutrons  can  be  used  to  produce  extra  fuel,  as  in  the  fast  breeder  reactor,  or  to  transmute 
long  half-life  waste  to  less  troublesome  isotopes,  such  as  in  the  Super  Phenix  reactor  near 
Cadarache  in  France,  or  in  some  combination  of  these  two  purposes  [9].  An  overview  of 
nuclear  energy  development  is  shown  in  Figure  5: 
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Figure  5.  Overview  of  Reactor  Technology  Development  [10] 


C.  FAST  SPECTRUM  NUCLEAR  REACTOR  DESIGN  FEATURES 

An  example  of  a  fourth  generation  system  currently  under  development  is  a  Lead- 
Cooled  Fast  Reactor  (LFR),  which  features  a  fast-spectrum  lead  or  lead/bismuth  eutectic 
liquid  metal-cooled  reactor  with  a  closed  fuel  cycle.  The  Pb  coolant  is  a  poor  absorber  of 
fast  neutrons  and  this  enables  the  realization  of  improved  sustainability  and  fuel  cycle 
benefits.  Pb  does  not  interact  vigorously  with  air,  water/steam,  or  carbon  dioxide, 
eliminating  concerns  about  exothermic  reactions.  It  has  a  high  boiling  temperature,  so  the 
prospect  of  boiling  or  flashing  of  the  ambient  pressure  coolant  is  realistically  eliminated. 
The  LFR  is  cooled  by  natural  convection  with  a  reactor  outlet  coolant  temperature  of  550 
°C,  possibly  ranging  up  to  800  °C  with  advanced  materials.  The  higher  temperature 
enables  the  production  of  hydrogen  by  thennochemical  processes  [11] 
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Figure  6.  Lead-Cooled  Fast  Reactor  Diagram  [10] 


Additionally,  due  to  its  neutron  economy,  an  LFR  can  be  designed  in  such  a 
manner  as  to  behave  like  a  converter  or  even  a  breeder  reactor,  therefore  enabling  the 
conversion  of  the  uranium  238  isotope  into  fissionable  plutonium  239.  A  neutron 
captured  by  uranium  238  results  in  the  formation  of  uranium  239,  which  has  a  half  life  of 
about  23  minutes  and  decays  into  neptunium  239  through  beta  decay.  Neptunium  239  has 
a  half  life  of  2.4  days  and  then  decays  into  plutonium  239,  also  through  beta  decay  [12]. 
The  added  benefit  of  the  high  neutron  flux  and  fast  spectrum  core  design  is  the  ability  of 
the  LFR  to  serve  in  the  actinide  management  role.  [13] 

Since  the  fission  of  fuel  materials  initially  present  in  the  core  creates  more  new 
fuel  (plutonium  239)  from  uranium  238,  as  a  side  product  of  fission,  a  properly  designed 
reactor  core  can  significantly  extend  the  time  between  refueling,  or  indeed  eliminate  it 
altogether.  An  advanced  reactor  design  that  maximizes  this  property  is  the  Small,  Sealed, 
Transportable,  Autonomous  Reactor  (SSTAR),  currently  under  development  by 
Lawrence  Livermore  and  Argonne  National  Laboratories  [11]. 
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D. 


SSTAR  PROGRAM  OVERVIEW 


The  objective  of  the  SSTAR  program  is  to  create  a  sealed  reactor  that  can  be 
delivered  to  a  site,  left  to  generate  power  for  up  to  30  years,  and  retrieved  when  its  fuel  is 
spent.  The  potential  for  nuclear  proliferation  would  be  minimized  by  sealing  the  cartridge 
core  inside  a  tamper  proof  cask.  The  reactor  would  be  monitored  and  operated 
autonomously  with  a  satellite  link  to  a  national  authority  or  international  authorities 
overseeing  the  reactors.  The  system  would  include  passive  safety  features  such  as  a  fast 
spectrum  core  with  a  strong  reactivity  feedback  that  enables  autonomous  load-following 
and  provide  passive  power  shutdown  in  case  of  a  loss-of-coolant  accident,  as  well  as 
natural  circulation  flow  of  the  Pb  coolant.  This  is  a  critical  issue,  as  the  reactor’s  primary 
market  is  thought  to  be  in  undeveloped  or  underdeveloped  regions  of  the  world,  where  it 
would  operate  with  almost  full  autonomy.  This  marketing  strategy  also  provides  the 
reasoning  for  the  small  size  of  the  reactor,  since  conventional  nuclear  stations  typically 
produce  about  a  gigawatt  of  electricity,  which  would  overwhelm  the  distribution  grid  in  a 
developing  country,  therefore  wasting  much  of  the  installed  power  [14]. 


TAMPER¬ 
PROOF  CASK 


STEAM 

GENERATOR 


Figure  7.  Diagram  of  a  Potential  SSTAR  Design  [14] 
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The  nominal  SSTAR  design  (20  MWe/  45  MWt)  is  18  meters  tall  and  3.2  m  in 
diameter,  which  enables  it  to  be  easily  moved  by  barge  or  rail.  The  producer  would  be 
responsible  for  delivering  the  sealed  unit  by  ship  and  truck  and  installing  it  at  the 
operating  location.  At  the  end  of  its  regular  operating  life,  the  old  reactor  would  be 
removed  for  recycling  or  disposal,  and  a  replacement  system  would  be  deployed.  [14] 


Figure  8.  SSTAR  Deployment  [14] 

SSTAR  would  be  coupled  with  an  S-CCfi  gas  turbine  Brayton  cycle  power 
converter,  which  enables  the  reactor  to  operate  with  a  higher  efficiency  when  compared 
to  the  traditional  Rankine  saturated  steam  cycle.  This  modification  is  enabled  by  the 
higher  core  temperatures  in  the  LFR  design,  attainable  with  the  primary  Pb  coolant  and 
highly  enriched  transuranic  (TRU)  nitride  fuel  in  15N  in  a  compact  core.  The  temperatures 
in  consideration  are  650°C  peak  cladding  temperature  and  561°C  core  outlet  temperature 
for  a  405  °C  core  inlet  temperature  [13].  However,  the  advantages  of  this  design  face  a 
critical  limitation  on  the  lifetime  of  this  and  other  Generation  IV  reactor  designs  from  the 
material  corrosion  which  results  from  radiation  damage  to  structural  materials 
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III.  MOTIVATION 


A.  RADIATION  DAMAGE 

Fission-emitted  neutrons  carry  with  them  up  to  5  MeV  of  energy.  Continuous 
neutron  bombardment  of  structural  materials  exposed  to  such  high  energy  neutron 
radiation  results  in  introduction  of  defects  within  the  material’s  crystalline  lattice.  The 
inelastic  collisions  between  the  impacting  neutrons  and  the  structural  atoms  cause  energy 
to  be  transferred  to  the  atoms.  Since  the  structural  materials  are,  for  the  most  part, 
comprised  of  stable  elements,  the  system  response  to  the  higher  energy  state  of  the  atom 
is  its  displacement  from  its  original  site  within  the  lattice.  This  disorder  within  the  perfect 
crystal  lattice  is  referred  to  as  a  Frenkel  disorder  [15].  It  results  in  the  formation  of  a 
Frenkel  pair  of  defects,  where  there  is  a  single  vacancy  and  a  single  interstitial  within  the 
structure  [15]. 


Figure  9.  Vacancy  Site  in  a  Crystalline  Structure  [16] 


Vacancies  are  sites  which  would  normally  be  occupied  by  an  atom  but  which  as  a 
result  of  the  displacement  are  unoccupied.  If  a  neighboring  atom  moves  to  occupy  the 
vacant  site,  the  vacancy  moves  in  the  opposite  direction  to  the  site  which  used  to  be 
occupied  by  the  moving  atom.  The  stability  of  the  surrounding  crystal  structure 
guarantees  that  the  neighboring  atoms  will  not  simply  collapse  around  the  vacancy.  In 
some  materials,  neighboring  atoms  actually  move  away  from  a  vacancy,  because  they  can 
better  fonn  bonds  with  atoms  in  the  other  directions  [15]. 
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Interstitials  are  atoms  which  occupy  a  site  in  the  crystal  structure  at  which  there  is 
usually  not  an  atom.  They  are  generally  high  energy  configurations.  There  are  two  basic 
kinds  of  interstitials:  Intrinsic  and  extrinsic  interstitials.  Intrinsic  interstitials  are 
interstitial  atoms  of  the  same  kind  as  the  atoms  of  the  crystal  ("self-interstitials").  They 
are  practically  nonexistent  in  elemental  crystals  (that  are  found  in  all  metals)  with  the 
significant  exception  of  Si,  where  intrinsic  interstitials  play  an  important  role  in  diffusion 
and  microdefect  formation.  Extrinsic  interstitials  are  interstitial  atoms  of  a  foreign 
(extrinsic)  type,  e.g.  C  in  Fe  or  O  in  Si  (“mixed  interstitials”).  They  may  diffuse 
directly  through  the  lattice  (i.e.,  without  the  help  of  vacancies)  and  play  an  important  role 
in  many  technically  relevant  materials  [15] 

o 


inters  rtf.]  ci( 
atom 

o 

Figure  10.  Interstitial  Atom  in  a  Crystal  Lattice  [16] 

Complexes  can  form  between  different  kinds  of  point  defects.  For  example,  if  a 
vacancy  encounters  an  impurity,  the  two  may  bind  together  if  the  impurity  is  too  large  for 
the  lattice.  Interstitials  can  fonn  'split  interstitial’  or  'dumbbell'  structures  where  two 
atoms  effectively  share  an  atomic  site,  resulting  in  neither  atom  actually  occupying  the 
site  [15].  These  interstitial  pairs  can  be  oriented  in  three  possible  directions: 

•  <100>  :  Atoms  displaced  in  only  in  one  axis  as  shown  in  Figure  11; 


14 


Figure  11.  <100>  Interstitial  [17] 

•  <1 10>  :  Atoms  displaced  in  two  axes  as  shown  in  Figure  12;  and 


Figure  12.  <1 10>  Interstitial  [17] 

•  <11 1>:  Atoms  displaced  in  all  three  axes  as  shown  in  Figure  3-8. 

<1 1 1> 


Figure  13.  <11 1>  Interstitial  [17] 


The  result  of  a  single  neutron  collision  is  not  a  single  defect;  on  the  contrary  the 
interaction  of  a  single  1  MeV  neutron  with  an  atomic  nucleus  transfers  up  to  about  100 
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keV  to  the  primary  knock-on  atom  (PKA)  of  iron.  Some  of  this  recoil  energy  is  lost  to 
interactions  with  the  electron  cloud,  resulting  in  a  somewhat  lower  kinetic  energy  that  is 
dissipated  in  atomic  collisions.  The  PKA  kinetic  energy  is  then  transferred  by  numerous 
subsequent  collisions  and  resultant  displacements,  producing  further  generations  of 
recoiling  atoms  at  lower  energies.  The  process  tenninates  when  the  kinetic  energy  of  the 
nth-generation  of  recoils  falls  below  that  needed  to  cause  additional  displacements  [18]. 
Closely  spaced  interstitials  and  vacancies  quickly  recombine  and  only  about  one-third  of 
the  initial  displacements  survive.  Typically,  this  leaves  a  vacancy-rich  cascade  core, 
surrounded  by  a  shell  of  interstitials.  The  majority  of  interstitials  quickly  cluster  to  form 
small,  disc-shaped  features  that  are  identical  to  small  dislocation  loops  [19].  Along  with 
interstitials,  these  loops  are  very  mobile.  Diffusion  of  individual  interstitials  and  loops 
within  the  cascade  region  causes  additional  recombination  prior  to  their  rapid  long-range 
migration  (unless  they  are  strongly  trapped  by  other  defects  or  solutes).  Although  they  are 
less  mobile  than  interstitials,  vacancies  also  eventually  diffuse.  Through  a  series  of  local 
jumps,  the  vacancies  and  solutes  in  the  cascade  quickly  begin  to  evolve  to  lower  energy 
configurations,  forming  small,  three-dimensional  clusters,  while  others  leave  the  cascade 
region  [19].  The  small  clusters  are  unstable  and  can  dissolve  by  vacancy  emission. 
However,  the  small  clusters  also  rapidly  diffuse  and  coalesce  with  each  other,  forming 
larger  nanovoids,  which  persist  for  much  longer  times.  Solute  atoms  bind  to  the  vacancies 
and  segregate  to  clusters.  The  vacancy  emission  rate  is  lower  from  vacancy-solute  cluster 
complexes.  Small  solute  clusters  remain  after  all  the  vacancy  clusters  have  finally 
dissolved.  Expressing  damage  exposure,  or  neutron  dose,  in  terms  of  displacements-per- 
atom  (dpa)  partially  accounts  for  the  effect  of  the  neutron  energy  spectrum  on  the 
generation  of  cascade  defects  and  the  net  residual  defect  production  scales  with  dpa  [20], 


B.  MACRO-LEVEL  EFFECTS  OF  RADIATION  DAMAGE 

This  time  evolution  of  microscopic  defects  results  in  the  following  emergent 
macroscopic  material  behaviors: 

•  Volumetric  Swelling:  irradiated  material  will  display  volumetric 
expansion  in  all  axes.  This  results  from  lattice  parameter  elongation  and 
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void  formation  [21].  An  example  of  radiation  induced  swelling  of 
autensitic  steel  is  shown  in  Figure  14 


Figure  14.  Radiation  Induced  Swelling  [21] 


•  Irradiation  Creep:  a  time-dependent,  constant  rate  mechanical 
deformation  of  material  occurring  slowly  at  stresses  below  the  ultimate 
tensile  stress.  Vacancies  cause  biased  absorption  of  interstitials  at 
dislocations  leading  to  net  dislocation  climb  [22].  An  example  of 
irradiation  creep  is  shown  in  Figure  15. 


Figure  15:  Irradiation  Creep  of  Nuclear  Fuel  Pin  Cladding  [23] 
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•  Radiation  hardening  and  embrittlement.  This  is  a  second  order 

phenomenon  which  is  caused  by  the  slipping  over  the  dislocation-linked 
clusters  of  point  defects.  This  leads  to  decreased  plasticity  of  the  irradiated 
material  and  ultimately  greater  proclivity  towards  rapid  catastrophic 
fracture  generation.  These  effects  become  even  more  pronounced  as  the 
ambient  temperature  increases  [18]. 

stress  concentration 


Figure  16:  Diagram  of  Material  Hardening  and  Embrittlement  Evolution  [18] 

C.  MATERIAL  ISSUES  FOR  GENERATION  IV  NUCLEAR  ENERGY 
SYSTEMS 

The  fact  that  nuclear  materials  are  subject  to  radiation  damage  is  well  known  from 
the  many  reactor  years  of  operation  accumulated  in  operation  of  Generations  I  and  II 
systems.  However,  when  we  compare  the  operating  conditions  present  in  early 
generations  of  nuclear  reactors  to  those  projected  for  Generation  IV  and  other  advanced 
nuclear  applications,  we  see  that  the  operating  environment  adversity  significantly 
increases. 
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Fission 

(Gen.  I) 

Fission 

(Gen.  IV) 

Fusion 

(Demo) 

NASA 
space  react. 

Structural  alloy 
maximum  temperaUue 

<300  C 

500-1  OOO’C 

550-1000T 

~1000°C 

Max  dose  for  core 
internal  structures 

-1  dpa 

-30-100  dpa 

-150  dpa 

-10  dpa 

Max  transmutation 
helium  concentration 

-0.1  appm 

~3- 1 0  appm 

-1500  appm 
(~  1 0000  appm 
for  SiC) 

-1  appm 

Coolants 

H,0 

He,  H,0.  Pb- 
Bi.  Na 

He.  Pb-Li,  Li 

Li,  Na,  or 
He-Xe 

Structural  Materials 

Zircaloy, 

stainless 

steel 

Ferritic  steel. 
SS, 

superalloys, 
C-  composite 

Ferritic/ 
martensitic 
steel,  V  alloy, 
SiC  composite 

Nb-IZr.  Ta 
alloys.  Mo 
alloys, 
superalloys 

Table  1.  Operating  Conditions  for  Advanced  Nuclear  Applications  [24] 


This  increase  in  adversity,  as  seen  from  Table  1,  make  the  materials  previously 
used  for  nuclear  applications,  namely  austenitic  steels  (such  as  316-L),  no  longer  suitable 
for  advanced  applications.  Therefore,  ferritic-martensitic  steels,  such  as  HT-9,  T-91  and 
EP823  have  to  be  considered.  The  chemical  composition  of  the  ferritic-martensitic  steels 
is  compared  to  an  example  of  autensitic  steel  in  Table  2. 


I.imii  il 
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Si 

Mu 

S 

P 

w 

Cr 

Ni 

Mo 

V 

Nl> 

N 

Ti 

Fe 

HT-9 

0.2 

0.25 

0.5 

.... 

.... 

0.5 

12 

0.56 

1 

0.3 

— 

— 

— 

84.7 

T-91 

0.1 

0.4 

0.45 

.... 

— 

— 

9 

— 

1 

0.2 

0.08 

0.05 

— 

88.7 

EPS23 

0.18 

1.05 

0.6 

— 

0.012 

0.65 

11.4 

0.7 

0.67 

0.4 

0.2 

— 

0.03 

84.1 

316-L 

0.035 

0.08 

2.0 

0.03 

0.040 

— 

16 

10 

2 

— 

— 

— 

— 

69.6 

Table  2.  Chemical  Composition  of  Potential  Generation  IV  Structural  Materials  [25] 


We  can  see  that  a  common  thread  among  the  materials  considered  is  that  they  are 
all  variants  of  iron-chromium  (Fe-Cr)  alloy,  but  the  following  are  the  practical  reasons 
why  these  materials  were  selected: 
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1 .  Unlike  stainless  steel  which  starts  to  swell  almost  immediately  after  being 
irradiated,  Fe-Cr  alloys  exhibit  a  delay  in  swelling,  wherein  they  do  not 
exhibit  that  property  until  receiving  100  dpa  of  radiation  damage. 

2.  Fe-Cr  alloys  perform  well  at  high  temperatures  and  pressures  required  for 
operation  of  advanced  nuclear  reactors,  and  they  have  favorable  rates  of 
corrosion  when  exposed  to  the  elements  intended  for  use  as  coolants  (Pb, 
Bi,  Na). 

3.  Their  welding  properties  are  well  known  and  there  is  experience  with  high 
temperature  helium  embrittlement  [26]. 

However,  despite  all  these  favorable  properties,  the  behavior  of  defect  clustering 
resulting  in  irradiation  creep,  hardening  and  embrittlement  is  poorly  understood  in  Fe-Cr 
alloys.  Considering  that  these  behaviors  are  known  to  occur  in  materials  used  in  earlier 
reactor  designs,  which  operate  in  much  more  benign  conditions,  it  becomes  imperative  to 
thoroughly  understand  them.  Additionally,  when  we  take  into  account  the  intended  life¬ 
time  of  a  Generation  IV  system  (around  60  reactor  years)  and  that  some  designs  (such  as 
SSTAR)  are  intended  to  operate  autonomously  for  extended  periods  of  time,  it  is  clear 
that  our  ability  to  model  the  time  evolution  of  Fe-Cr  alloys  is  critical  for  safe  operation  of 
these  systems. 

D.  MULTI-SCALE  APPROACH  TO  MATERIAL  SCIENCE  MODELING 

Since  radiation  induced  defects  occur  at  the  atomistic  level,  the  most  accurate 
method  of  calculating  their  energetic  properties  is  to  solve  the  Schrodinger  equation  (or 
Dirac  equation,  if  relativistic  effects  are  important)  in  its  many-body  electron  fonn. 
However,  this  is  an  exceptionally  difficult  and  computationally  intensive  process,  as  the 
currently  “standard”  model  for  condensed  matter  physics,  the  Density  Functional  Theory 
(DFT)  using  Local  Density  Approximation  (LDA),  is  currently  limited  to  100-1000 
atoms.  The  largest  DFT  simulation  to  date  is  1080  B  atoms  (with  3840  electrons 
simulated)  on  Lawrence  Livermore  National  Laboratory’s  2000  CPU  Linux  cluster  [23]. 
While  the  DFT  approach  is  generally  successful  in  predicting  structures  and  macroscopic 
properties,  it  under-predicts  band  gap  energies,  over-predicts  lattice  parameters,  and 
predicts  wrong  ground  states  for  some  magnetic  systems  (e.g.,  Fe).  Nonetheless,  for  the 
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purposes  of  material  and  reactor  engineering,  requiring  time  evolution  of  macro  scale 
behaviors  such  as  void  swelling,  hardening,  embrittlement,  creep,  stress  corrosion 
cracking,  the  first-principles  calculation  is  not  suitable  [27].  Therefore,  in  order  to  arrive 
at  a  successful  time-evolved  simulation  of  a  finite  structural  element,  several  layers  of 
approximation  are  necessary. 

The  multi-scale  modeling  approach  to  material  science  is  only  natural,  because 
the  introduction  of  defects,  their  evolution  and  mutual  interaction  and  ultimately  their 
effect  on  the  mechanical  properties  of  the  material  all  occur  on  very  broad  timescales,  as 
can  be  seen  in  Figure  17. 
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Figure  17.  Temporal  and  Spatial  Placement  of  Phenomena  [23] 


The  modeling  of  these  phenomena  in  a  global  sense  is  a  chain  of  events  in  which  the 
results  of  each  previous  step  serve  as  the  source  of  inputs  and  basis  of  comparison  for  the 
subsequent  steps.  Fitting  the  model  at  the  next  higher  level  of  approximation  is  done  by 
comparing  its  results  for  properties  best  modeled  by  the  method  on  a  finer  scale.  By  using 
a  bottom-up  approach,  based  on  first  principles  and  built  upon  scale  reversible  models 
(when  possible)  starting  from  the  sub-nanometer  scale  (where  the  building  blocks  of 
matter  are  established,  hence  providing  material  unity  and  technology  integration), 
complete  characterization  of  properties  in  materials  and  processes  at  different 
dimensional  and  time  scales  can  be  achieved.  In  addition,  a  sequence  of  experiments  is 
designed  to  follow  the  modeling  effort  and  validates  the  simulation  result.  A  visual 
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Timescale 


description  of  the  overall  multi-scale  modeling  effort,  complete  with  supporting 
experiments  is  shown  in  Figure  18. 
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Figure  18.  Global  View  of  Multi-scale  Material  Modeling  Effort  [28] 
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IV.  THEORY  AND  METHODOLOGY  FOR  COMPUTATION 


A.  INTRODUCTION 

Modem  materials  science  evolved  from  two  basic  elements— metallurgy  and 
alchemy— and  has  over  the  period  of  the  last  couple  of  decades  greatly  expanded  its  field 
of  study.  Its  natural  complexity  had,  for  a  very  long  time,  kept  its  progress  in  check  as 
the  difficult  physical  concepts  could  only  be  applied  to  idealized,  pure  forms. 
Nonetheless,  Moore’s  Law  compounding  of  computer  power  has  recently  enabled 
numerical  simulation  of  increasingly  complex  systems.  While  the  most  exact  calculations 
are  based  on  first  principles  of  quantum  mechanics,  alternatively  referred  to  as  ab  initio 
calculations,  their  sheer  complexity  has  yet  to  be  conquered  by  Moore’s  Law.  The 
complexity  of  these  calculations  is  clear  to  anyone  who  has  ever  taken  an  introductory 
course  in  quantum  mechanics,  which  culminates  in  solving  the  Schrodinger  Equation  for 
a  hydrogen  atom.  However,  since  most  useful  materials  contains  significantly  more 
complex  elements,  the  difficulty  is  compounded  as  it  becomes  necessary  to  solve 
Schrodinger’s  many-electron  equation  in  the  following  form: 

n*  =  [t  +  v  +  u\  a-  =  +  £>«)  + !>«,?,)  *  =  £# 

-  4  i  J  (1) 

Where  H  is  the  electronic  molecular  Hamiltonian,  N  is  the  number  of  electrons  and  U  is 
the  electron-electron  interaction.  The  operators  T  and  U  are  so-called  universal  operators 
as  they  are  the  same  for  any  system,  while  V  is  system  dependent  or  non-universal.  In 
essence,  the  greatest  complexity  in  ab  initio  methods  arises  from  the  many-body 
interaction  term  U.  While  there  are  many  clever  methods  for  reducing  the  computational 
power  required  to  solve  this  equation  for  a  complex  system,  those  remain  a  discussion  for 
another  day.  A  powerful  approach  has  been  developed  in  Ref.  [29],  that  uses  classic 
interatomic  potentials  and  incorporates  ab  initio  and  experimental  results,  in  order  to 
enable  the  study  of  crystalline  defects  and  their  long  range  interactions  on  a  much  larger 
scale. 
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B. 


THE  MANY-BODY  POTENTIAL 


The  models  generally  used  for  this  work  are  colloquially  known  as  “many-body” 
potentials,  and  can  be  globally  grouped  in  two  categories:  the  embedded-atom  models 
and  the  second  moment  approximation  [30].  However,  the  issue  facing  material  modelers 
for  advanced  nuclear  applications  is  that  the  overwhelming  body  of  work  in  this  field  is 
focused  on  either  pure  elements  or  intennetallic  compounds;  only  a  few  authors  address 
concentrated  alloys. 

The  methodology  used  for  this  work  was  developed  by  Dr.  Alfredo  Caro,  et  al. 
of  the  Lawrence  Livermore  National  Laboratory,  first  published  in  the  Physical  Review 
Letters  article  of  August  2005  [29].  The  objective  of  this  methodology  is  to  address 
arbitrarily  complex  systems  of  concentrated  alloys  with  complex  heat  of  formation.  What 
makes  this  work  especially  appropriate  is  that  this  methodology  is  applied  to  Fe-Cr  alloys 
which  are  of  interest  for  advanced  nuclear  applications. 

The  “many-body”  potentials  are  based  on  the  summation  of  atom  energies 
which  provides  the  total  energy  of  the  system.  The  atom  energies  are  themselves 
composed  of  two  contributions;  namely,  embedding  and  pair  potential  tenns.  For 
heteroatomic  systems,  such  as  binary  alloys  composed  of  elements  A  and  B,  this  reads: 


where  a  and  P  stand  for  elements  A  and  B  sitting  at  sites  i  and  j  within  the  crystalline 
lattice,  F’s  are  the  embedding  functions  for  either  type  of  element,  and  V’s  and  p’s  are  the 
pair  potentials  and  densities  between  a-P  pairs.  The  functions  pap  and  Vap  therefore 
describe  the  properties  of  the  alloy.  The  variety  in  expressions  of  embedding  energies, 
densities,  and  pair  potentials  encompasses  a  great  similitude  of  models. 

The  greatest  issue  in  application  of  this  model  to  the  case  of  Fe-Cr  is  that  unlike 
other  binary  alloys,  its  fonnation  energy  is  not  a  symmetric  function.  Instead,  in  the  case 
of  Fe-Cr  this  function  is  highly  nonsymmetrical  and  it  even  changes  sign  at  low  Cr 
composition  [35],  as  seen  in  figure  19. 


24 


Figure  19.  Ab  Initio  Calculations  of  Fe-Cr  Mixing  Enthalpy  [35] 

Therefore,  in  contrast  to  the  binary  alloys  with  symmetric  formation  energy 
such  as  Au-Ni  and  Fe-Cu,  the  standard  approach  using  a  cross  pair  potential  term  is  not 
good  enough  to  reproduce  the  properties  of  the  systems.  By  using  the  potentials  already 
described  in  the  literature,  adjusting  the  alloy  tenn  in  equation  (2)  and  focusing  on 
nonlinearities  built  upon  the  pair  potential  cross  term  alone,  this  methodology  manages  to 
create  a  model  which  successfully  departs  from  ideality  and  matches  the  heat  of 
formation  behavior  as  described  by  Par  Olsson  et  al.  in  [35]  and  as  shown  in  Figure  20. 


Figure  20.  Formation  Energy  of  the  Alloy  as  Predicted  by  the  Potential  Used  in  This 

Work  [29] 

C.  THE  FREE  ENERGY  EXPRESSION 

The  free  energy  model  uses  an  effective  representation  of  two  pure  element 
potentials  as  reported  by  Mendelev  for  Fe  [36]  and  Wallenius  for  Cr  [37],  with 
normalized  densities,  which  for  a  =  A,  B  reads 
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Pa  Pa!  £a,eq> 

=  K(e°J  -  K'(e 

VaJr)  =  V'Jr)  +  2F"'(e^)p"(r), 

W/ 

where  the  superscript  0  stands  for  original,  f  °a  for  the  density  on  a  lattice  site  at 


equilibrium  MPa Xr,y9)]’  anc*  the  prime  ’  for  derivative.  These  transformations  do  not 


alter  the  properties  of  the  pure  elements  but  have  the  advantage  of  minimizing  the 
contribution  of  the  embedding  tenn  to  the  formation  energy  of  the  alloy,  thus  enabling 
unrelated  pure  element  potentials  to  be  combined  in  this  alloy  description,  as  shown  by 
the  free  energy  of  a  random  solution  alloy  with  composition  x  at  temperature  T  in 
equation  4: 

g(x,  T)  =  gKf(x,  T)  +  gmixU  T)  +  AgU,  T), 

(4) 

where  gref  is  the  compositional  weighted  free  energy  of  the  pure  components,  gmix  is  the 
free  energy  contribution  from  the  entropy  of  mixing  for  a  random  alloy,  and  Ag  is  the 
excess  Gibbs  energy  of  mixing,  which,  when  expressed  by  a  Redlich-Kister  expansion 
[38]  reads 


A g(x,  T)  =  jc(1  -  jc)  Lp(T)(  1  -  2 x)p, 

p= o 


(5) 


where  Lp  is  the  p-th  order  binary  interaction  parameter.  This  parameter  is  also  the  object 
of  two  major  simplifications  in  this  methodology,  as  it  is  generally  a  function  of 
temperature.  Therefore,  in  order  to  reduce  the  complexity  of  the  description  we  neglect 
the  excess  vibrational  entropy  and  assume  that  the  formation  energy  does  not  depend  on 
T.  This  leads  to  a  simplified  version  of  equation  (5). 


Ag(jt,  T)  =  A //(a)  =  a(1  -  x)  ]T  Lp(  1  -  2  x)p. 

p= o 


(6) 


The  factors  for  the  Redlich-Kister  expansion  come  from  the  Par  Olsson  ab  initio 
calculations  [35]  and  are  given  in  the  table  below. 
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L0 

Li 

l2 

L3 

U 

0.41566 

0.0814134 

-0.0101899 

0.267659 

-0.248269 

Table  3. Values  of  the  Red 

lich-Kister  expansion  coefficients  for  Eq.  (6) 

D.  THE  CROSS  PAIR  POTENTIAL 


The  functional  form  of  the  cross  potential  is  based  on  the  analytic  mode  of  the 
alloy  in  which  the  species  that  sits  on  site  i  can  be  either  A  or  B,  but  both  are  embedded 
in  the  same  average  environment,  as  discussed  by  Ackland  and  Vitek  [31]: 


£rand  =  4  X  VM(r,j)  +  ,v|  X  VBB(ru) 

+  2xAxB  £  VAB(r,j)  +  xAFA(p)  +  xBFB(p) 


(7) 


where  the  p  is  defined  as 

P  = 

Therefore,  the  contribution  of  the  embedding  tenns  to  the  mixing  energy,  AEemb’  is 

=  xA[ Fa(p)  -  FA(p  =  I )] 

+  xB[FB(p)  -  FB(p  =  1)J 


Using  a  Taylor  expansion  of  F  around  p=  1  and  using  equation  (3),  it  becomes  clear  that 


this  contribution  is  quadratic  in  {p  -  1),  and  therefore  small  for  small  variations  in  p- 1 
as  seen  in 

A£emb  =  ,\AFA(p  =  I  Hp~  1) 2+XBF"(p  =  1  )(p  ~  I)2.  (10) 

The  specific  potentials  used  for  this  study  reduce  the  embedding  term  contribution 
to  the  fonnation  energy  down  to  ~1  meV/atom  at  x  =  0.5,  making  it  negligible  in 
comparison  to  the  target  value  of  -100  meV  for  Fe-Cr  alloys  [35].  This  enables  us  to 
write  down  the  energy  contribution  from  the  pair  potential  terms  as  follows, 

A£Pair  =  a(1  -  x){2vab  -  {vA  +  vB)}. 
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It  is  therefore  assumed  that:  a)  the  nonlinearity  is  dependent  only  on  the  pair 
potential  and  b)  Vab  is  a  function  of  both  x  and  r,  separable  as  a  product  of  h(x)uAB(r), 
with  the  form, 

VabU.  r)  =  Mx^V^r)  +  VM(r)].  ([2) 

The  h(x)  is  also  a  4th  order  Reidlich-Kister  polynomial,  whose  factors  are  obtained 
through  global  minimization  of  equation  (7)  using  Mathemathica™  ,  and  listed  in  Table 
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Ho 

hi 

h2 

h3 

h4 

0.883644 

-0.059302 

0.644634 

-1.342524 

0.918932 

Table  4.  4th  Order  Polynomial  Coefficient  for  h(x),  as  Obtained  by  Minimization 


This  methodological  approach  yields  a  formation  energy  curve  indistinguishable 
from  the  target  function  obtained  by  the  ab  initio  calculations  of  Par  Olsson  [35],  as  seen 
in  Figure  19.  Therefore,  it  can  be  concluded  that  it  is  reasonable  to  utilize  this  potential 
for  modeling  the  formation  energies  of  point  defects  in  Fe-Cr  alloys. 


E.  THE  MOLECULAR  DYNAMICS  CODE 


According  to  Ohno  et  al.  [42],  the  method  of  classical  molecular  dynamics  (MD) 
was  first  proposed  by  Alder  and  Wainwright  in  1957.  The  essence  of  the  approach  comes 
from  two  assumptions:  a)  that  there  is  a  set  of  defined  interatomic  force  functions  or  a 
potential,  as  we  have  shown  in  the  beginning  part  of  this  chapter;  and  b)  the  numerical 
ensemble  in  which  the  number  of  particles  N,  the  temperature  T  and  the  chemical 
potential  p  are  constant.  If  the  above  two  statements  are  correct,  the  equation  of  motion 
for  the  atoms  is  the  usual  Newtonian  equation 


n2r  n 

m  '  =  F  F  =  -V  V  V 

1  dt2  1,1  'tT  y 


(13) 


where  m;  is  the  mass  of  the  i-th  particle,  r,  its  coordinate  and  Fi  the  force  acting  on  it.  The 
code  solves  this  equation  numerically  and  performs  energy  minimization  of  the  ensemble, 
therefore  achieving  thennal  equilibrium  after  a  sufficient  number  of  iterations.  This  in 
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turn  provides  us  with  the  energetics  of  a  sample  in  which  all  the  atoms  have  achieved 
their  least  energetic  state  and  are  in  a  “relaxed”  state. 

One  of  the  greatest  limitations  of  the  MD  approach  has  been  in  coping  with  the 
effects  of  finite  system  size  and  presence  of  surfaces.  In  order  to  reduce  these  effects,  the 
code  uses  periodic  boundary  conditions.  With  this  approach,  all  the  particles  are  placed 
inside  a  box,  or  a  unit  cell.  If  the  particle  goes  outside  the  cell,  it  is  brought  back  in  from 
the  opposite  side  of  the  cell.  While  this  approach  works  marvelously  with  perfect  crystal 
systems,  introducing  a  disturbance  such  as  a  point  defect  creates  problems  because  the 
code  is  inclined  to  see  a  “mirror  image”  of  the  defect  across  the  periodic  boundary 
conditions,  therefore  unintentionally  creating  a  cascading  effect  of  mutual  influence. 

One  technique  for  minimizing  the  effects  of  finiteness  of  the  simulated  system 
is  to  sum  the  forces  exerted  on  the  i-th  particle  from  all  other  particles  inside  a  certain 
“sphere  of  influence”  defined  by  the  cut-off  radius  Rc  centered  on  the  particle.  Of  course 
this  summation  must  be  performed  even  if  the  particles  are  not  in  the  same  unit  cell,  but 
in  the  image  cells  [43].  Ideally,  having  an  adequately  large  sample  size  in  which  the 
sphere  of  influence  radius  around  the  i-th  atom  does  not  encompass  any  atoms  affected 
by  the  presence  of  its  mirror  image  solves  this  problem. 
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V.  PROCEDURE  AND  RESULTS 


A.  RESEARCH  PLAN 

The  thesis  research  was  conducted  according  to  the  following  schedule: 

•  Task  1:  The  Study  of  Defects  in  Pure  Elements.  The  objective  of  this  section  of 
the  thesis  is  to  prove  that  the  Caro  potential  [29]  adequately  approximates  results 
obtained  by  ab  initio  calculations.  This  will  be  done  by  calculating  formation 
energies  of  point  defects  in  pure  element  samples  and  comparing  those  results 
against  values  available  in  literature,  obtained  through  ab  initio  calculations. 

•  Task  2:  Evaluating  the  Vacancy  Energy  of  Formation.  This  section  is  aimed  at 
evaluating  the  behavior  of  vacancy  formation  energy  in  samples  with  increasing 
chromium  concentration.  This  will  serve  to  define  the  specific  energetics  of  this 
point  defect  for  a  range  of  Cr  concentration  in  relation  to  a  linear  interpolation 
between  vacancy  fonnation  energies  in  pure  elements. 

•  Task  3:  Obtaining  the  Formation  Energies  of  Interstitials.  In  this  section  the 
research  will  focus  on  defining  the  behavior  of  formation  energies  as  a  function  of 
Cr  concentration  for  a  variety  of  interstitial  configurations,  as  well  as  exploring 
their  relative  relationships. 

Put  together,  the  scope  of  these  results  will  serve  to  prove  the  applicability  of  this 
type  of  a  semi-empirical  potential  for  approximating  the  results  of  ab  initio 
calculations  and  provide  static  and  dynamic  calculations  on  a  scale  currently 
unobtainable  by  ab  initio  simulation. 
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B.  CHARACTERIZING  THE  PERFORMANCE  OF  THE  CODE:  THE 

STUDY  OF  DEFECTS  IN  PURE  ELEMENTS 

In  order  to  benchmark  the  performance  of  the  newly  developed  potential  [29] 
against  other  available  Fe  and  Cr  potentials  and  methodologies,  it  was  necessary  to 
compare  its  prediction  of  formation  energies  for  the  12  point  defects  types  of  interest: 

•  Single  Vacancy  in  a  pure  Fe  sample 

•  Single  Vacancy  in  a  pure  Cr  sample 

•  Self-interstitial  (Fe-Fe)  <1 10>  in  a  pure  Cr  sample 

•  Self-interstitial  (Fe-Fe)  <1 1 1>  in  a  pure  Fe  sample 

•  Self-interstitial  (Fe-Fe)  <1 10>  in  a  pure  Fe  sample 

•  Self-interstitial  (Fe-Fe)  <  1 1 1>  in  a  pure  Cr  sample 

•  Mixed  interstitial  (Fe-Cr)  <1 10>  in  a  pure  Fe  sample 

•  Mixed  interstitial  (Fe-Cr)  <1 1 1>  in  a  pure  Fe  sample 

•  Self-interstitial  (Cr-Cr)  <1 10>  in  a  pure  Cr  sample 

•  Self-interstitial  (Cr-Cr)  <  1 1 1>  in  a  pure  Cr  sample 

•  Mixed  interstitial  (Fe-Cr)  <1 10>  in  a  pure  Cr  sample 

•  Mixed  interstitial  (Fe-Cr)  <1 1 1>  in  a  pure  Cr  sample 

The  benchmarking  could  only  be  performed  at  the  endpoints  of  the  Cr  concentration 
diagram,  since  no  other  physically  sound  models  for  the  Fe-Cr  alloys  could  be  found  in 
literature. 

The  simulation  samples  were  prepared  by  manually  modifying  the  1024  atom, 
perfect  crystal  samples  of  pure  Fe  and  pure  Cr.  For  the  vacancy  cases,  this  modification 
consisted  of  deleting  the  atom  located  at  coordinates  (0,0,0)  in  both  samples.  For  the 
interstitial  cases,  these  modifications  consisted  of  adding  an  additional  atom  (atom  1025) 
to  the  end  of  the  sample  and  shifting  the  first  atom  on  the  list  (atom  1).  The  location  of 
atom  1025  was  determined  by  adding  0.5  A  to  the  coordinates  in  the  axes  where  the 
displacement  was  to  occur,  i.e.,  in  x  and  y  axes  for  <11 0>  type  interstitials  and  in  all  three 
axes  for  <11 1>  type  interstitials.  The  value  of  0.5  A  was  subtracted  from  the  respective 
coordinates  of  atom  1,  therefore  creating  the  relevant  interstitial  pair.  The  value  of  0.5  A 
was  detennined  by  trial  and  error,  as  initial  attempts  to  use  a  larger  initial  interstitial 
spacing  (0.9  A)  produced  results  in  which  the  interstitial  pair  would  fly  apart  during 
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energy  minimization.  For  self-interstitials,  an  atom  of  the  same  type  would  be  inserted 
into  the  sample,  while  the  mixed  interstitial  samples  received  an  atom  of  the  opposite 
type. 

The  samples  were  then  run,  using  LLNL’s  Multiprogramatic  Capability  Cluster  - 
MCR.  This  is  a  large  (1 1.2  TF)  tightly  coupled  Linux  cluster  with  1,152  nodes,  each  with 
two  2.4-GHz  Pentium  4  Xeon  processors  and  4  GB  of  memory.  MCR  runs  the  LLNL 
CHAOS  software  environment,  which  incorporates  the  Red  Hat  Linux  operating  system 
[44].  The  Molecular  Dynamics  (MD)  code,  Moldy  [45],  was  set  up  for  constant 
temperature  of  0.01  K,  no  pressure,  no  chemical  potential,  with  free  sample  volume,  and 
a  run  of  3000  steps  with  a  time  step  of  10'5  seconds.  This  number  of  steps  was  necessary 
in  order  to  ensure  that  the  sample  had  ample  time  to  relax  around  the  induced  defect  and 
that  the  final  energy  values  converged.  At  the  end  of  the  run,  the  final  total  enthalpy 
values  were  extracted  for  each  output  and  compared  against  those  obtained  from  running 
the  perfect  crystal  samples  under  the  same  conditions.  The  table  listing  the  final  values  of 
total  enthalpy  for  each  sample  is  shown  in  Table  5. 


Sample 

Total  Enthalpy  [eV] 

1024  Fe  perfect  crystal 

-4221.373163 

1 024  Cr  perfect  crystal 

-3928.353757 

1023  Fe  vacancy 

-4215.461833 

1 023  Cr  vacancy 

-3921.989661 

1025  self-interstitial  <11 0>  (Fe-Fe)  in  Iron 

-4221.974874 

1025  self-interstitial  <11 0>  (Fe-Fe)  in  Chromium 

-3929.176430 

1025  mixed  interstitial  <1 10>  (Fe-Cr)  in  Iron 

-4222.064124 

1025  mixed  interstitial  <11 0>  (Fe-Cr)  in  Chromium 

-3928.756038 

1025  self-interstitial  <1 1 1>  (Fe-Fe)  in  Iron 

-4221.546602 

1025  self-interstitial  <1 1 1>  (Fe-Fe)  in  Chromium 

-3929.176430 

1025  mixed  interstitial  <1 1 1>  (Fe-Cr)  in  Iron 

-4221.896441 

1025  mixed  interstitial  <1 1 1>  (Fe-Cr)  in  Chromium 

-3928.756244 

1025  self-interstitial  <11 0>  (Cr-Cr)  in  Chromium 

-3926.664855 

1025  self-interstitial  <1 1 1>  (Cr-Cr)  in  Chromium 

-3926.579530 

Table  5.  Table  of  Final  Values  of  Total  Enthalpy  for  Pure  Element  Samples 
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The  formation  energies  for  the  vacancy  (Efy  )  in  pure  elements  are  then  computed 
from  the  total  enthalpy  of  the  pure  element  sample  ( Eff  )  and  the  total  enthalpy  of  the 

sample  with  the  vacancy  ( Eff  ),  according  to  the  following  fonnula: 


-*1024 


1023 

1024 


(14) 


The  fonnation  energies  for  self  interstitials  (Eg)  were  computed  using  a  very 
similar  approach,  but  incorporating  the  fact  that  samples  which  contain  the  interstitial 
consist  of  1025  instead  of  1023  atoms,  and  therefore  correlating  the  energies  in  the 
following  way: 

77*  pure 

EFI  =  EmZZ1  -1025-(^gp,  (15) 


where  f 


int  errstitial 
1 025  atoms 


is  the  total  enthalpy  of  the  sample  with  1025  atoms  and  a  self  interstitial 


and  EZMoms  is  the  total  enthalpy  of  the  sample  with  the  pure  element  in  pure  crystal 


1 024  atom  arrangement. 

For  the  mixed  interstitial  samples,  however,  a  different  methodology  must  be 
applied,  defining  the  formation  energy  of  the  interstitial  as  follows: 


EFI  - 


E 


ICr  int  erstital 
1025  atoms 


-E 


reference 


where 


E  =  E 

reference  ^  int  erpolation 


+  E 


from—hof 


and 


(16) 


E  t  ,,  =  nr  •  Epuref  +nF  ■ EpureFf 

interpolation  Cr  per  atom  be  per  atom 


Efrom-hof  >  however,  is  a  product  of  a  polynomial  fit  to  Ah[eV/atom]  as  seen  in 

figure  3-4,  and  as  we  are  working  in  the  region  below  6  %  of  Cr  concentration,  where 
this  difference  function  is  negative,  we  must  only  use  the  Ah  values  from  a  series  of 
samples  with  an  increasing  Cr  concentration,  relaxed  by  running  them  through  the  Moldy 
code,  as  seen  below: 
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xCr 

Ah[eV/atom] 

0 

0 

0.00097 

-9.65E-05 

0.006836 

-0.000564056 

0.012695 

-0.000492174 

0.015625 

-0.000884431 

0.019531 

-0.001196044 

0.030273 

-0.000627474 

0.038086 

-0.000667051 

0.047852 

1.16076E-05 

Table  6.  Values  of  Ah[eV/atom]  as  Obtained  by  Modeling 


When  fitted  with  a  4th  order  polynomial  using  Excel,  this  produces  the  following  heat  of 
formation  function  in  the  sub-6%  Cr  concentration  region: 

_ _  E-Eideal  [mev/atl 

0.0002  1  J 


-0.0012  • 


-0.0014 

Figure  21.  Heat  of  Formation  Fit  Below  6%  Cr  concentration 


Therefore,  when  x  is  calculated  as  the  value  of  x  =  1/1025  and  the  corresponding  y  turns 
out  to  be  -7.85209  X  10'05  meV/atom.  E,  ,  .  is  the  1025th  multiple  of  that  value  of  y 

from-noj  r 

and  is  then  entered  into  (16). 

The  results  of  this  characterization  are  displayed  in  Table  7. 
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Table  7.  Formation  Energies  for  the  Pure  Element  Characterization 


C.  EVALUATING  THE  VACANCY  ENERGY  OF  FORMATION 

In  order  to  analyze  the  effect  of  increased  Cr  concentration  on  the  formation 
energy  of  a  single  vacancy,  it  was  necessary  to  characterize  this  behavior  with  a  series  of 
samples  of  increasing  Cr  concentration.  Also,  since  it  is  unknown  what  the  influence  of 
the  distance  between  the  solute  Cr  atoms  and  the  vacancy  is,  it  was  necessary  to  analyze 
the  entire  ensemble  of  possible  configurations. 

Therefore,  11  samples  of  1024  atoms  with  Cr  concentrations  ranging  from  a 
singular  Cr  atom  to  17  %  were  simulated  using  the  molecular  dynamics  code  Moldy  [45], 
and  set  up  for  a  constant  temperature  of  0.01  K,  no  pressure,  no  chemical  potential,  with 
free  sample  volume,  and  a  run  of  3000  steps  with  a  time  step  of  10'5  seconds.  These  runs 
provided  the  reference  values  of  total  enthalpy  for  the  pure  crystal  samples.  Those  values 
are  shown  in  Table  8. 
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Cr 

Concentration 

Total  Enthalpy 

[eV] 

9.77E-04 

-4221.45017 

0.01563 

-4221.88415 

0.03027 

-4221.0655 

0.04785 

-4219.7705 

0.06055 

-4217.80019 

0.07422 

-421418631 

0.08984 

-4210.24463 

0.09668 

-4209.23258 

0.10645 

-4206.00442 

0.11523 

-4203.40938 

0.17188 

.4184.4663 

Table  8.  Total  Enthalpy  of  the  reference  samples 


Thereafter,  each  of  the  samples  was  run  through  a  sequence  of  two  programs 
using  a  LINUX  script.  The  first  program  went  through  the  sample  and  removed  one  atom 
at  a  time,  creating  a  vacancy.  Only  Fe  atoms  were  removed  in  order  to  maintain  the  Cr 
composition.  At  the  end  of  operation,  approximately  1000  vacancy  samples  were  ready, 
representing  all  possible  configurations.  Next,  all  of  the  newly  created  samples  were 
simulated  using  the  molecular  dynamics  code  Moldy,  and  set  up  for  a  constant 
temperature  of  0.01  K,  no  pressure,  no  chemical  potential,  with  free  sample  volume,  and 
a  run  of  3000  steps  with  a  time  step  of  10'5  seconds.  In  the  end,  the  LINUX  code 
extracted  the  final  values  of  total  enthalpy  from  each  sample  output  file. 

A  histogram  of  energies  could  then  be  produced  for  each  set  of  samples 
containing  the  same  number  of  Cr  atoms.  Over  this  histogram,  a  Gaussian  distribution 
function  was  fitted  using  ORIGIN™  with  the  canonical  form  of 


y  =  To  + 


x—xc 

w 


(17) 
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where  yo,  A,  xc  and  w  are  fitting  parameters.  In  our  case,  xc  provides  us  with  the  mean 
value  of  the  total  enthalpy  of  the  sample  and  w  provides  us  with  the  width  of  the  2/3 
width  of  the  Gaussian,  which  represents  one  standard  deviation  of  results  and  gives  us  the 
width  of  the  error  bars.  Each  mean  value  for  total  enthalpy  of  a  sample  including  a 
vacancy  was  then  used  in  conjunction  with  the  respective  value  of  total  enthalpy  of  the 
perfect  crystal  to  provide  the  fonnation  energy  of  the  vacancy  at  the  given  Cr 
concentration.  Since  this  equation  is  first-order  linear,  the  uncertainty  in  ,  as 

obtained  from  the  Gaussian  fit  (w),  transfers  directly  into  Efy.  The  vacancy  formation 
energies  obtained  are  shown  in  Table  9. 


Cr 

Concentration 

Vacancy  Energy  of 

Formation  (eV) 

w/2 

9.77E-04 

-1.713903803 

0.00215 

0.01563 

-1.715881926 

0.0749 

0.03027 

.1.71444999 

0.01249 

0.04785 

-1.722512512 

0.01765 

0.06055 

-1.761504721 

0.05799 

0.07422 

-421414418.6 

0.06145 

0.08984 

-1.785316452 

0.06531 

0.09668 

-1.799338475 

0.06517 

0.10645 

-1.838959668 

0.06851 

0.11523 

-1.812067625 

0.06319 

0.17188 

-1.863116579 

0.05324 

Table  9.  Vacancy  Formation  Energy  as  a  Function  of  Cr  Concentration 

D.  OBTAINING  THE  FORMATION  ENERGIES  OF  INTERSTITIALS 

The  procedure  for  obtaining  the  formation  energies  of  self-interstitials  is  quite 
analogous  to  those  described  for  the  evolution  of  vacancy  formation  energy  and  for  self¬ 
interstitial  cases  in  the  pure  elements.  After  obtaining  the  reference  values  for  the  samples 
in  the  perfect  crystal  state,  the  samples  are  operated  on  by  a  program  which  places  an 
interstitial  atom  added  at  each  and  every  atom  site  containing  a  Fe  atom.  The  program 
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displaces  the  original  atom  by  0.5  A  away  from  its  original  location  and  adds  its  mirror 
image  atom.  The  displacement  takes  place  in  x  and  y  axes  for  interstitials  <1 10>  and  all 
three  axes  for  interstitials  <1 1 1>.  After  creating  approximately  1000  samples  for  each  Cr 
concentration  examined,  the  script  then  runs  each  of  them  in  a  Moldy  code,  with  the 
inputs  set  at  a  constant  temperature  of  0.01  K,  no  pressure,  no  chemical  potential,  with 
free  sample  volume,  and  a  run  of  3000  steps  with  a  time  step  of  10'5  seconds.  The 
extracted  values  of  final  enthalpy  are  then  sorted  in  an  energy  histogram  and  a  Gaussian 
fit  is  applied,  thus  providing  us  with  the  mean  energy  value  for  each  selected 
concentration  as  well  as  an  indication  of  uncertainty  based  on  the  configurational 
characteristics  of  the  sample.  The  values  of  final  enthalpy  for  pure  crystal  samples, 
median  values  for  those  containing  an  interstitial,  their  individual  uncertainties  and  the 
self-interstitial  fonnation  energy  calculated  via  equation  (15)  are  shown  in  Table  10. 


Cr 

Formation  Energy 

w/2 

Formation  Energy 

w/2 

Concentration 

<110>  eV 

<110> 

<11 1>  eV 

<111> 

0.000976 

3.519643423 

0.004155 

3.942243423 

0.00729 

0.006829 

3.516942656 

0.006155 

3.941462656 

0.00883 

0.012683 

3.512976597 

0.019815 

3.932626597 

0.014195 

0.030244 

3.506226279 

0.031625 

3.911206279 

0.020885 

0.047805 

3.504111963 

0.04858 

3.844751963 

0.07725 

0.074146 

3.494965576 

0.106945 

3.811195576 

0.06379 

0.106341 

3.498952173 

0.154455 

- 

- 

0.12 

3.48620728 

0.171505 

- 

- 

0.145366 

3.47881083 

0.19146 

- 

- 

0.171707 

3.47969771 

0.19536 

- 

- 

Table  10.  Energetics  of  Self-Interstitial  Samples 


To  obtain  the  evolution  of  mixed  interstitial  fonnation  energies,  we  apply  an 
almost  exact  procedure,  except  that  when  inserting  an  atom  we  are  introducing  a  Cr, 
rather  than  a  Fe.  As  this  changes  the  global  composition  of  the  sample,  we  must  correct 
for  this  using  equations  (16).  In  order  to  obtain  Efrom  hof  value  it  is  necessary  to  make  a 
dual  zone  polynomial  fit.  In  this  case  zone  A  is  below  6%  Cr  concentration  and  zone  B 
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The  resulting  energetics  of  mixed  interstitials  are  shown  below  in  Table  1 1 : 


Cr 

Formation  Energy 

Formation  Energy 

Concentration 

<110>  eV 

w/2  <110> 

<111>  eV 

w/2  <11 1> 

0.001951 

3.428325907 

0.007 

3.596965907 

0.004155 

0.007805 

3.433025476 

0.01334 

3.591255476 

0.021995 

0.013659 

3.432935366 

0.0165 

3.603105366 

0.033825 

0.03122 

3.467658628 

0.06466 

3.608058628 

0.062045 

0.04878 

3.509552851 

0.0963 

3.628022851 

0.12394 

0.075122 

3.550068635 

0.105995 

3.633258635 

0.15496 

0.107317 

3.546490688 

0.10822 

3.621460688 

0.184035 

0.120976 

3.537463852 

0.109665 

3.604573852 

0.192875 

0.146341 

3.530354706 

0.126305 

3.546284706 

0.140005 

0.172683 

3.484627439 

0.124385 

3.501687439 

0.18874 

Table  11.  Mixed  Interstitial  Energetics  in  eV 
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VI.  ANALYSIS  AND  CONCLUSION 


A.  VACANCIES 

As  we  can  see  from  Figure  24,  our  results  for  vacancy  formation  energy  in  Fe  are 
lower  than  DFT  calculations..  However,  when  we  consider  the  fact  that  most  of  these 
calculations  were  either  thermally  unequilibrated  (unrelaxed)  or  had  a  volume  constraint, 
this  is  becomes  a  much  more  believable  result.  Both  of  these  conditions  led  to  additional 
“frustration”  within  the  sample,  and  accordingly  to  a  higher  fonnational  enthalpy. 
References  for  all  values  quoted  are  listed  at  the  end  of  the  chapter. 

On  the  other  hand,  when  we  compare  the  vacancy  formation  energy  in  Cr  with 
other  results  (Figure  25),  as  well  as  those  calculations  which  produce  both  results  (Figure 
26),  we  see  that  it  falls  neatly  in  between  ab  initio  and  other  classical  results.  From  both 
of  these  comparisons  we  can  conclude  that  our  potential  gives  a  reasonable  prediction  of 
vacancy  fonnation  energies  for  both  pure  elements,  as  well  as  for  their  mutual 
relationship. 

The  more  extensive  analysis  of  vacancy  formation  energy  behavior  as  a  function 
of  Cr  concentration  yields  some  fascinating  results.  Prior  to  this  work  it  has  been 
supposed  that  the  vacancy  formation  energy  would  follow  a  direct  interpolation  between 
the  pure  element  values.  As  we  can  see  in  Figure  27,  the  true  formation  energy  deviates 
from  the  assumed  linear  interpolation  in  Figure  26  in  the  region  below  6%  Cr 
concentration,  which  is  the  same  region  in  which  the  heat  of  formation  curve  shows  the 
same  deviant  behavior  from  linear  interpolation.  In  real  terms,  this  finding  means  that  it  is 
more  likely  for  vacancies  to  form  in  alloys  with  a  Cr  concentration  smaller  than  6%,  then 
those  above  this  value.  The  exact  nature  of  this  effect  on  the  physical  properties  of  the 
alloy  is  as  of  now  undeterminable;  however,  it  can  be  hypothesized  that  it  will  be 
observable.  Additionally,  the  increasing  spread  of  values  as  a  function  of  Cr 
concentration  shows  a  strong  interactivity  between  the  impurities  and  vacancies. 
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Figure  24.  Vacancy  Fonnation  Energies  in  Iron 


This  Work  WaUenius  -  EAM,  Dudarev,  DFT,  Olsson,  DFT,  Const  Experiment 
Free  V,  Relaxed  Const.  V,  Unrelaxed  V,  Unrelaxed 


Figure  25.  Vacancy  Formation  Energies  in  Chromium 
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2.7 


2.64 

2.56 


2.1 


Figure  26.  Linear  Interpolation  Between  EFV  in  Fe  and  Cr 


Figure  27.  Evolution  of  Vacancy  Formation  Energy  as  a  Function  of  Cr 
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B. 


SELF-INTERSTITIALS 


As  we  can  see  from  Figures  28  and  29  our  results  match  both  the  relative  order 
and  the  magnitude  of  target  ab  initio  calculations.  The  critical  part  of  this  fitting  was 
achieving  the  same  relative  order  of  values  for  mixed  and  self  interstitials  in  iron,  as  both 
theory  and  experiment  suggest  the  mixed  interstitial  to  be  the  prime  defect  in  Fe-Cr  alloys 
[45]. 

Examining  the  plot  of  interstitial  formation  energy  as  a  function  of  Cr 
concentration  yields  some  exceptionally  interesting  infonnation.  While  the  formation 
energy  of  the  interstitial  in  <1 10>  remains  relatively  stable  around  3.5  eV  per  defect,  and 
actually  drops  lower  as  Cr  concentration  increases,  the  <1 1 1>  curve  shows  some  highly 
unstable  behavior.  An  initial  examination  of  energy  histograms,  produced  by  compiling 
and  sorting  the  final  output  enthalpies  of  all  samples  initialized  with  a  <1 1 1>  interstitial, 
showed  a  behavior  that  deviated  significantly  from  that  which  was  expected.  For 
example,  Figure  30  shows  the  energy  histogram  for  samples  containing  10.6  %  Cr  atoms 
and  an  interstitial  in  <110>.  The  energies  in  the  histogram  follow  a  nonnal  distribution 
which  can  be  approximated  and  analyzed  with  a  Gaussian  curve  fit  as  overlaid.  For 
samples  which  contain  just  a  single  Cr  atom  the  values  are  also  normally  distributed  (see 
Figure  31);  however,  their  variance  is  much  smaller  because  most  configurations  interact 
equally  with  the  interstitial.  The  only  deviants  are  those  where  the  Cr  atom  is  in  the  near 
neighborhood  of  the  interstitial  site. 

However,  when  we  take  a  look  to  the  histogram  of  sample  enthalpies  with  0.006 
Chromium  concentration,  in  Figure  32,  we  make  the  astonishing  observation  that  there  is 
not  one,  but  two  distinct  distributions  of  energies.  Upon  closer  examination,  we  notice 
that  the  larger  distribution  has  a  mean  value  which  is  equivalent  to  that  of  the  distribution 
for  samples  with  the  <1 10>  interstitial  at  the  same  Cr  concentration.  This  inevitably  leads 
us  to  conclude  that  a  number  of  samples,  even  though  initially  seeded  with  a  <11 1> 
interstitial,  converted  to  an  interstitial  in  the  <11 0>  direction  during  the  process  of  energy 
minimization.  Physically,  this  means  that  the  <1 1 1>  interstitial  is  not  a  preferred  state  for 
the  defect  and  even  if  induced  it  will  tend  to  convert  into  a  <110>  type  interstitial.  The 
fraction  of  samples  which  convert  into  a  <11 0>  type  interstitial  increases  steadily  as 
shown  in  Figure  33.  This  increase,  combined  with  the  fact  that  the  differential  between 
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self-interstitial  formation  energies  in  <11 0>  and  <11 1>  decreases  as  seen  in  Figure  34, 
makes  it  impossible  to  discern  the  existence  of  a  peak  with  final  enthalpies  for  samples 
containing  the  <11 1>  type  interstitial  after  approximately  10  %  because  the  variance  has 
increased  so  much  that  the  two  distributions  completely  overlap. 

C.  MIXED  INTERSTITIALS 

When  looking  at  the  mixed  interstitials  and  their  formation  energy  evolution  as  a 
function  of  Cr  concentration,  we  note  several  important  differences  when  compared  to 
the  self-interstitials. 

In  this  case,  it  is  the  <11 0>  interstitial  that  has  a  lower  energy  of  formation, 
therefore  being  the  preferable  state.  However,  both  mixed  interstitial  curves,  as  can  be 
inferred  from  comparing  Figure  35  with  Figure  34,  are  above  those  for  self-interstitials. 
However,  they  are  both  exceptionally  stable  and  well  behaved,  so  much  so  that  they  run 
almost  parallel  above  10  %  in  concentration  and  even  have  very  little  overlap  when  their 
variances  are  taken  into  account. 
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Formation  Energy  Mixed  Formation  Energy  Mixed  Formation  Energy  Self  Formation  Energy  Self 
Interstitial  <  1 1 0>  Interstitial  <  1 1 1  >  Interstitial  <  1 1 0>  Interstitial  <  1 1 1  > 


Figure  28.  Formation  Energies  for  Interstitials  in  Iron  [45] 
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Cr  self-interstitial  <11 1> 


Figure  29.  Formation  Energies  for  Self-Interstitials  in  Chromium 


Data:  Count5_Count 
Model:  Gauss 

Equation:  y=yO  +  (A/(ui*sqrt(P  l/2)))'exp(-2,‘((x-xc)AD )"2) 
Weighting: 
y  No  weighting 
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Figure  30.  Energy  Histogram  for  Self-Interstitial  Samples  in  the  <11 0>  direction  at 

higher  Cr  Concentration 
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Figure  31.  Histogram  of  Energies  for  Self-Interstitial  Samples  in  the  <11 0>  direction  at 

lower  Cr  Concentration 
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Figure  32.  Histogram  of  Energies  for  Self-Interstitial  Samples  in  <1 1 1> 
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Figure  33.  Evolution  Self-Interstitial  Formation  Energies 
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Figure  35.  Evolution  of  Mixed-Interstitial  Formation  Energies 
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Figure  36.  Combined  Interstitial  Evolution  Plot 
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D.  CONCLUSIONS 

If  we  look  at  Figure  36,  we  can  see  that  all  interstitial  types  behave  in  a  very 
similar  manner  higher  the  Cr  concentration.  It  is  not  clear  what  physical  significance  this 
prediction  could  have,  but  it  would  have  to  be  validated  through  either  experimental 
methods  or  ab  initio  calculations,  before  further  discussion  can  ensue.  Nonetheless,  the 
ability  of  EAM  to  offer  such  thorough  evaluations  should  be  noted. 

To  conclude,  it  has  been  proven  that  this  methodology  [29]  adequately 
approximates  results  obtained  through  ab  initio  calculations.  It  has  also  been 
demonstrated  that  by  using  this  methodology,  the  studies  of  point  defects  can  be  of  a 
larger  scale  and  more  exhaustive  than  is  the  case  with  the  currently  available  ab  initio 
models. 

The  results  produced  in  this  research  will  go  into  studies  of  dynamic  properties  of 
point  defects  such  as  vacancy  migration  and  defect  clustering.  In  the  meantime,  though,  a 
interim  step  is  a  study  of  precipitation  properties  of  Fe-Cr  alloys,  which  will  have  a  high 
fidelity  considering  that  the  underlying  assumptions,  i.e.,  the  interatomic  potentials,  have 
already  been  validated  as  providing  an  accurate  approximation  of  ab  initio  results,  and  by 
inference  real  world  physical  behavior. 

In  order  to  study  this  behavior,  a  hybrid  Monte  Carlo-molecular  dynamics  code 
named  MCCASK  was  developed  by  A.  Caro  and  B.  Sadigh  in  2005  at  Lawrence 
Livermore  National  Laboratory.  The  code  performs  sequences  of  Monte  Carlo  events  and 
Molecular  Dynamics  time  steps.  In  this  way,  the  equilibrium  concentrations  in  the  alloy 
are  obtained. 

The  molecular  dynamics  part  of  the  MCCASK  code  is  based  in  MDCASK,  a 
molecular  dynamics  code  that  was  first  developed  to  study  radiation  damage  in  metals. 
MDCASK  solves  the  equation  of  motion  of  the  atoms  in  the  sample;  energy  and  forces 
on  each  atom  are  calculated  using  an  interatomic  potential,  and  the  equations  of  motion 
are  integrated  to  obtain  the  next  values  of  positions  and  velocities.  The  Monte  Carlo  part 
of  the  new  code  MCCASK  corresponds  to  a  parallel  MC  code  in  the  transmutation 
ensemble  (T,P,N,Dm),  with  N  the  total  number  of  atoms  and  Dm  the  difference  in 
chemical  potential  [46]. 
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Computational  materials  science  is  a  young  field,  which  has  just  recently  been 
handed  the  tools  of  its  trade— high-end  computing  hardware  and  software— and  over  the 
coming  years  it  is  inevitable  that  major  advances  in  modeling  will  result  in  real  life 
impact,  not  only  in  the  field  of  nuclear  energy,  but  in  everyday  applications.  The  future  is 
bright. 
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VII.  SUMMARY  OF  RESEARCH 


In  summary,  it  can  be  said  that  the  suitability  of  Fe-Cr  alloys  for  use  in  cladding, 
wrappers  and  in-core  structural  components  in  Generation  IV  advanced  nuclear  systems 
necessitates  a  thorough  understanding  of  their  behavior  under  irradiation,  to  include 
formation  energies  of  point  defects  such  as  vacancies  and  interstitials  in  the  steel  matrix. 
This  work  used  a  powerful  methodology  developed  by  A.  Caro,  et  al.  [29]  to  determine 
the  energetics  of  vacancies  and  interstitials  both  in  pure  elements  and  as  a  function  of  Cr 
concentration.  An  overview  of  the  methodology  is  shown  in  Figure  37. 
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Figure  37.  Research  Methodology  [29] 
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Results  of  this  research  proved  that  the  methodology  is  appropriate  and  can  be 
utilized  for  large  scale  simulation,  and  that  it  approximates  ab  initio  results  satisfactorily. 
In  addition  to  that,  the  evaluation  of  point  defect  formation  energy  as  a  function  of  Cr 
concentration,  gave  some  interesting  predictions  concerning  their  static  properties  and 
stability.  A  finding  relating  the  stability  of  the  <11 1>  Fe  self-interstitial  to  local  Cr 
concentration  is  summarized  in  Figure  38. 


xCr 

Figure  38.  Development  of  the  Stability  Function  for  Fe  Self-interstitial  in  <1 10> 

We  can  therefore  say  that  the  results  obtained  in  this  thesis  can  be  used  for  studies 
of  dynamic  properties  and  time  evolution  of  point  defects.  The  stability  predictions  will 
need  to  be  verified  through  ab  initio  or  experimental  data,  but  at  the  very  least  they  can 
serve  to  demonstrate  the  potential  benefits  of  using  molecular  dynamics  and  semi- 
empirical  potentials  to  approximate  first  principles  calculations. 
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